An M E B D F Package for the Numerical Solution of Large Sparse Systems of Stiff Initial Value Problems

K e y w o r d s M E B D F , Large sparse systems of stiff IVPs and differential-algebraic equations, Time dependent PDEs. 1 . I N T R O D U C T I O N In th i s pape r , we will be concerned wi th the numer ica l so lu t ion of the in i t ia l value p rob l e m dy M~x = f (x ,y) , y(xo) = Yo. (1.1) We will be p a r t i c u l a r l y concerned wi th two i m p o r t a n t cases. (i) M = I and so we have a sys t em of in i t ia l value p rob lems which we will a s sume to be stiff. (ii) T h e m a t r i x M is s ingular so t h a t we have a sys tem of d i f ferent ia l -a lgebra ic equa t ions (DAEs) . Most , if no t all, of t he efficient numer ica l me thods for the solut ion of stiff sys t ems of t he form (1.1) a re impl ic i t and , as a d i rec t resul t , i t is necessary to solve a nonl inear sys t em of a lgebra ic equa t ions a t each t ime s t ep in o rder to compu te t he requi red numer ica l solut ion. In t he p resen t pape r , we will be concerned wi th the case where t he Jacob ian m a t r i x ~ is ve ry large and sparse. Such The authors are very grateful to W. Schiesser and his coworkers for their continuing interest in this project and for making their codes available to us. We benefitted greatly from free access to their library subroutines for space discretization. We are also grateful to two referees for their many helpful comments which greatly improved this paper. 0898-1221/01/$ see front matter (~) 2001 Elsevier Science Ltd. All rights reserved. Typeset by .AA/~S-TEX PII: S0898-1221 (01)00137-7 122 T.J. ABDULLA et al. systems can arise quite naturally in practical applications. They also arise from the use of the method of lines (MOL) in the solution of time dependent PDEs. If the PDE being solved has just one space dimension, then the system of ODEs to be solved is normally banded and this is relatively easy to solve using standard software (an efficient banded MEBDF [1] solver is already available from the web page of one of the present authors [2]). If, however, the PDE has more than one space dimension, then the resulting system of ODEs is generally sparse and this presents a much more formidable problem to initial value solvers. A crucial component in developing an efficient sparse ODE solver is, of course, the derivation of an effective solver for the large sparse system of nonlinear algebraic equations occurring at each time step. There are at least two different approaches to this problem. The first is to use an iterative technique such as Gauss-Seidel, SOR, or Krylov methods [3]. This approach is very efficient as regards storage space, can often be implemented efficiently on a parallel computer, particularly if Jacobi iteration is used, and may often be the only realistic option in three space dimensions. However, the obvious disadvantage is that in some cases, it may experience convergence difficulties. The second approach is to base the algebraic equation solver on a Newton iteration scheme and to use a standard sparse equation solver to solve the resulting system of linear algebraic equations. Such an approach will normally require more storage than is required by an iterative solver, due to the possibility of fill in, but will enjoy the well-known advantages of a direct, as opposed to an iterative, solver. A major disadvantage of this approach is that there may be convergence difficulties with the Newton iteration. In this paper, we will concentrate on the use of a direct solver. However, our code is written in a modular form so that it is relatively straightforward for a user to replace the direct solver by an iterative one if so required. Most direct solvers for sparse linear algebraic equations are based on Gaussian elimination and they usually employ a pivoting strategy which is a compromise between pivoting solely for stability and pivoting to avoid fill in. In this paper, we will concentrate on two sparse solvers, namely the Yale solver YSMP [4] and the Harwell solver MA28 [5]. As we shall see, we choose the Yale solver for ODEs and the Harwell solver for differential-algebraic equations (DAEs). Another very efficient solver is the special purpose solver developed in [6]. We have not implemented it in our code because it requires the Jacobian to be stored as a full N x N matrix and this is not practical for the problems we are interested in. In the next section, we will briefly describe the MEBDF approach for differential-algebraic equations. A detailed analysis can be found in [7]. 2. M O D I F I E D E X T E N D E D B A C K W A R D D I F F E R E N T I A T I O N F O R M U L A E Modified extended backward differentiation formulae (MEBDF) [1] were originally proposed for the numerical solution of stiff initial value problems (IVPs) in an at tempt to derive a class of multistep integration formulae which have better stability properties and higher orders of accuracy than standard backward differentiation formulae (BDF). The basic MEBDF algorithm was subsequently modified to produce the code MEBDFDAE [2] which is applicable to linearly implicit differentiai-algebraic equations of index < 3, as well as to stiff IVPs. In an extensive comparison of codes [8], it was found that MEBDFDAE compares very favourably with certain other state of the art codes on a large set of challenging problems. Further evidence of the good performance of the MEBDF codes can be found on the web page of one of the present authors [2]. At the present time, the main options allowed by the MEBDFDAE code are to specify either full or banded Jacobians which are computed either numerically or analytically. The purpose of the present paper is to extend these options to allow the solution of large sparse systems of initial value problems both for ODEs and DAEs. Having done this, our aim is to compare the performance of the sparse version of MEBDF with that of the sparse BDF code LSODES [9]. In particular, we will examine what the more stable, higher-order but more expensive MEBDF MEBDF Package 123 approach has to offer in the method of lines (MOL) solution of time dependent PDEs after the PDE has been semidiscretized in space. In what follows, we will explain the MEBDF algorithm for differential-algebraic equations. For ease of presentation, we will consider the one-step case, that is, the case where only one past value Yn is used. The general MEBDF algorithm consists of the following three stages. STAGE 1. Assuming the data (xn, Yn) are known, we compute Yn+l from M (Yn+l -Yn) = hS ( X n * l , Y n , l ) , using a Newton iteration scheme. Once convergence to a solution Yn+l has been obtained from the Newton iteration, we compute _, Yn+z Yn and Yn+I h S (Xn'-t-l,Yn÷l) ~-MYln÷l • STAGE 2. Compute Yn+2 from M (Yn+2 Yn+l) ---h S (Xn+2,Ynq-2)" Similarly to Stage 1, once convergence has been obtained to Yn+2' we compute -, Yn+2 Yn+I and Yn+2 = h f (xn+2, Yn+2) = Myln+2 • STAGE 3. Compute Yn+l from M Yn+l Yn -~ [ f (XnT2,Yn+2) -Jr f (XnTl,Yn+l)] -'~ h f (Xn+l, Yn+l). Note that this formula is the one-step MEBDF of order 2. This is explained, for example, in Halter and Wanner [1, pp. 267-270]. As a first approximation to Yn+l and f ( X n + l , y n + l ) , we use Yn+l and f ( x n + z , Y n + l ) from Stage 1. This then completes one step forward of the MEBDF algorithm. For reasons explained more fully in [10], Stages 1 and 3 are normally relatively cheap compared with Stage 2 and so the total work in implementing a full MEBDF step, consisting of the above three individual stages, is normally not much more than for a single BDF step. Indeed, if the Newton iteration in Stage 3 converges in just one iteration, as is often the case, then Stage 3 requires no extra function evaluations, but simply calls for the solution of a linear system